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Abstract. We first analyse the effect of a square root transformation to the time variable 
on the convergence of the Crank-Nicolson scheme when applied to the solution of the heat 
equation with Dirac delta function initial conditions. In the original variables, the scheme is 
known to diverge as the time step is reduced with the ratio, A, of the time step to space step 
held constant and the value of A controls how fast the divergence occurs. After introducing 
the square root of time variable we prove that the numerical scheme for the transformed 
partial differential equation now always converges and that A controls the order of conver- 
gence, quadratic convergence being achieved for A below a critical value. Numerical results 
indicate that the time change used with an appropriate value of A also results in quadratic 
convergence for the calculation of the price, delta and gamma for standard European and 
American options without the need for Rannacher start-up steps. 

Keywords. Heat equation; Crank-Nicolson scheme; convergence; Black-Scholes; European 
option; American option; asymptotics; time change. 



1. Introduction 



The Crank-Nicolson scheme is a popular time stepping scheme for the numerical approx- 
imation of diffusion equations and is particularly heavily used for applications in computa- 
tional finance. This is due to a combination of favourable properties: its simplicity and ease 
of implementation (especially in one dimension); second order accuracy in the timestep for 
sufficiently regular data; unconditional stability in an L2 sense. 

However, there are well-documented problems which ar ise from the fact that the scheme is 
at the cusp of stability in the sense that it is A-stable (see Smithl . [l985l ) but does not dampen 
so-called stiff components effectively (i.e., is not strongly A-stable). Specifically, a straight- 
forward Fourier analysis of the scheme applied to a standard finite difference discretisation of 
the heat equation shows that high wave number components are reflected in every time step 
but asymptotically do not diminish over time. This gives rise to problems for applications 
with non-smooth data, where the behaviour of these components over time plays a crucial 
role in the smoothing of solutions. Examples where this is releva,nt incl ude Dirac initial data, 
as they appear naturally for adjoint equations (see Giles &: Siili . 20021 ) . and piecewise linear 
and step function pavoff data in the computation of option prices fsee Wilmott. Howison & 
Dewynne, 19951 ). 

There, the situation is e xacerbated if sensitivities of the solution to its input data (so-called 
'Greeks') are needed; see IShaw (|l999f ) for an early discussion of issues with time stepping 
schemes in this context. 

There is a sizeable body of literature co ncerned with schemes with better stability prop- 
erties, in particular L- stable schemes (see Smith . 19851 ) which share with the underlying 
continuous equation the property that high wave number components decay increasingly fast 
in time. Examples for schemes with this property include sub-diagonal Fade schemes based 
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on rational approximations of the exponential function. In Rannacher ( 19741 ) . Rannacher 
proposes to replace the first 2jU steps of a higher-order Fade scheme of order (e.g., 
the Crank-Nicolson scheme for = 1) by a low-order sub-diagonal Fade scheme with order 
{li — (e.g., backward Euler for fj. = 1), and shows that this restores opt i mal o rder for 
diffusion equations; the c ase // = 1 is already analys e d in [Luskin &: Rannacher (1982). 

It is demonstrated in iFoolev. Vetzal k Forsvthl l20o£ how this procedure can be used 
to obtain stable and accurate approximations to option values and their sensitivities for 
various n on-smooth payoff functions. T his has been widely adopted in financial engineering 



practice. Khaliq. Voss &: Yousuj ( 2007 ). Khaliq extend this to more general Fade schemes 



and demonstrate numerically the stability and accuracy in practice for options with exotic 
payoffs. Further ad aptations are possible, for instance, where Wade. Khaliq, Yousuf, Vigo- 
Aguiar & Deininger (|2007l l give an application to discretely sampled barrier options, where 
disconti nuities are i ntrod uced at certain points in time and a new 'start-up' is needed (see 
already Rannach er^ ('l974j) for the provision of such restarts in an abstract context). 

Carter &: Giles (2007) analyse the behaviour of the Crank-Nicolson numerical scheme when 
used to solve a convection-diffusion equation with constant coefficients and Dirac delta func- 
tion initial conditions. The numerical solution obtained by this method diverges as the time 
step goes to zero with the ratio A = k/h of the time step, k, to the space step, h, held 
constant. Their analysis in the frequency domain shows that there are errors associated with 
high wave numbers which behave as 0{h~^) and which will eventually overwhelm the 0{h?) 
errors associated with low wave numbers as the time step goes to zero, thus explaining why 
the errors in this scheme will eventually increase as the time step is reduced. Keeping the 
ratio A constant is desirable because the Crank-Nicolson central difference scheme is second 
order consistent in both space and time and the numerical scheme is more efficient if this 
property is exploited. The authors extended their analysis to show why the incorporation 
of a small number of initial fuHv irnplicit time steps fi.e.. the Rannacher scheme. Luskin & 



Rannacher ( 19821 ): Rannacher! ( 19741 )) eliminates this divergence. 

In this paper, an alternative method of avoiding the divergence is proposed and analysed. 
The idea is to introduce a time change into the partial differential equation (FDE) by trans- 
forming the original time variable, t, to 

i = Vt 

and solving the equation numerically using the Crank-Nicolson scheme in the new variables. 
The time change will be applied to the heat equation ([1]) on M x [0, T] and can be considered 
natural as the heat equation with suitable initial data admits similarity solutions which depend 
on xl\/t. From a probabilistic angle, the heat equation is the Kolmogorov forward equation 
for the transition density of Brownian motion whose standard deviation at time t is 
The main result of the paper is given by the following Theorem. 

Theorem 1. The Crank-Nicolson central difference scheme with uniform time steps exhibits 
second order convergence (in the maximum norm) for the time-changed heat equation, with 
Dirac initial data, as the time step k tends to zero with X = k/h held fixed, provided that 
A < l/\/2- For A > l/\/2 the scheme is still convergent with order 1/A^. 

The peculiar dependence of the convergence rate on the mesh ratio is in fact seen to be 
sharp up to a logarithmic factor. It also follows from the analysis of the heat equation that 
the computational complexity for an optimal choice of the mesh ratio is lower than for the 
Rannacher scheme in its optimal refinement regime. 
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As an added benefit, we obtain experimental se cond order converg e nce fo r the value, delta 
and gamma of European and American options. Forsvth &: Vetzall ( 20021 ) observe that at- 



the-money prices of European and American options computed by a Crank-Nicolson finite 
difference scheme exhibit a reduced convergence order of 1 in the time step k, which improves 
to 2 for Rannacher start-up in the case of European options, but only to 3/2 for American 
options. This last observation is rationalised by the square-root behaviour of the value and 
exercise boundary for short time-to-expiry. A heuristic adaptive time-stepping scheme based 
on this observation is shown there to restore second order convergence. We show here, by 
means of numerical tests, that the square root time change introduced above provides a 
similar remedy even without Rannacher start-up. This is not surprising bearing in mind the 
relation between the time change and a non-uniform time mesh in the original time variable. 
Precisely, the time-changed scheme with constant steps of size VT/N, where N is the total 
number of time steps for time up to T, is equivalent to a non-uniform scheme in the original 
time variable with time points tm = {km)"^ for m = 0, \fT jk. The step size hence increases 
in proportion to 2\ft and the smallest (the first) step is of size T/N"^. 

The remainder of this article is organised as follows. In Section [21 we apply the time 
change transformation to the heat equation, describe the transformed numerical scheme, and 
calculate the discrete Fourier transform of the numerical solution. Then, in Section [31 we 
perform asymptotic analysis of the error between the transform and the transform of the true 
solution, identifying four wave number regimes which will then be used in the later analysis. In 
Sectional we use the results of the asymptotic analysis to determine the asymptotic behaviour 
of the errors between the numerical solution and the true solution, from which we deduce the 
convergence behaviour of the transformed scheme and hence Theorem [H Section [5] contains 
numerical results illustrating these findings and compares the computational complexity to 
the Rannacher scheme. In Section [6l the time change transformation is applied to the Black- 
Scholes equation and the solution method is used to calculate the gamma of a European call 
option. The behaviour of the resulting errors is described and explained in terms of the values 
of the strike and volatility. Finally, in Section [71 the time change transformation is applied 
to a pen altv method used to calculate the price of an American out, as described in Forsvth 
& Vetzal (|2002h . Section (81 concludes. 



2. The Crank-Nicolson scheme applied to the transformed heat equation and 

ITS behaviour in the Fourier domain 

To simplify the analysis, we restrict attention to the heat equation (i.e., we do not consider 
any convection terms) 

1 

Ut = -Uxx, (1) 

for t G [0, T], X G M, where the solution, n, satisfies ti(x,0) = 5(x), where 5 is the Dirac delta 
function and which, after the change of time variable, becomes 

11^ — t Uxx ) 

where i = \/t. 

The transformed equation is then discretised using the Crank-Nicolson method, leading to 
the following numerical scheme 

k h? +2* h? ' 
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where C/" is the numerical solution at the spatial grid node j at the 'time' step n, where the 
space grid, in principle, extends from — oo to +00, with Xj = jk for j G Z. In the transformed 
'time' direction, we have f = ik for 1 < i < where is the number of timesteps. We can 
simplify this to 

where A = k/h. 

To study the behaviour of this sch eme in the Four i er dom ain we apply the d iscrete Fourier 
transform to this equation (e.g., see lCarter &: GilesI (j2nn7l ^ and lStrand ll9H(h ). Multiplying 



equation ([2]) by e**^^ , summing over j and simplifying, we obtain the following recurrence for 
the Fourier transform 

j=+oo 
j=-oo 

at successive time steps, 

1 + 2(n + 1)A2 sin2 (^^^^ [7"+i(s) = (^1 - 2nA2sin2 (^^^^ 

We note that U^{s) = 1, as the initial condition is a delta function at (0,0), so we have 

frN(. ^ l(l-0(l-20...(l-(jV-l)g) 

(1+0(1 + 26. ..(1 + iVO ' 

where ^ = 2A^ sin^(s/i/2) and we want to study the behaviour of U^{s) as — )• 00 with kN 
and k/h held fixed. 

In the subsequent analysis, it is often simpler to describe a condition on s, by stating the 
related condition on ^, but it should be noted that the relationship between s and ^ depends 
on h. 

In what follows, and to simplify the notation, we will fix T = 1, so that 

iv = i = J-. 

k hX 

We note that the exact solution of the PDF is given by 

and its Fourier transform is given by u{s,t) = exp {-sH/2). So, for i = i = 1, we have 

u{s, 1) = exp (-s'^/2). 

As in ICarter fc Gilei (|2007l ^ we can estimate the errors in the numerical solution (at time 
t = i = 1), i.e. the differences between the values of Uj^ and u{xj, 1), by applying the inverse 

Fourier transform to C/^(s) — u{s,l). 

3. Asymptotic analysis of the Fourier transform 

In this section we identify four wave number regimes for the Fourier transform and ob- 
tain useful asymptotic estimates for C/^(s), in each of these regimes. Initially, we start by 
identif ving three regimes, the low, intermediate and high wave number regimes fas in Carter 
& Giles~5jr 

but then find it useful to subdivide the high wave number regime into two 



further regimes. 
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In essence, the wave number regimes are determined by the locations of the real zeros of 
the equation 

U'^is) = 0, 

which are at Sm, for m = m* , . . . , N — 1, for some integer m* > 1, which depends only on A, 
where 

n\2 ■ 2 ' Sm.h\ 1 

2 A sm 



2 y m 

and we note that m* > Note that each Sm is a function of h and — ^ oo as /i — )• 0, 
so these wave number boundaries increase in absolute terms while always being less than or 
equal to vr/Zi. 

The low wave number regime (which we will refer to as regime I) is, in part, defined by the 
requirement that all the terms in the numerator of the expression for U^{s) in equation 1^ 
should be positive. This requires that ^ < 1/{N — 1). If we assume the stronger condition 
that ^ < we can write a truncated Taylor expansion for 

N-l N 

log{U^{s)) = ^ log(l -mC)-Y, log(l + ^0 

m=l m,=l 

from wh ich we can later d erive an asymptotic estimate for U^{s). 



As in ICarter &: GilesI (|2007l ). the low wave number regime is further restricted by the 



condition that s < , for a value r such that r < ^, which ensures that the remainder terms 
in the Taylor expansion go to zero as /i — )• 0, so that we can derive a useful approximation for 
U^{s) — u{s, 1). As this condition also ensures that ^ < asymptotically, this condition 
alone is sufficient to define the low wave number regime. 

Next, there exists an intermediate wave number regime (regime II ) defined by the range of 
wave numbers for which ^ < but which are not in regime I (cf. ICarter k GilesI (|2007h ). 



The high wave number regime is defined by the requirement that ^ > If this condition 

on ^ is satisfied then some of the terms in the numerator of the expression for [/^(s) in 
equation ^ will be negative. There will be an integer m > such that ^ G [l/(m + 1), 1/m] 
and the number of positive terms will remain fixed as N increases while the number of 
negative terms will increase with A^. A particular negative factor (1 — r^) can be rewritten as 
— r^(l — with < 1, and we can then write a truncated Taylor expansion for C/^(s) 

in terms of 1/^ which will then lead to an asymptotic value for C/^(s). 

At first sight, it might seem necessary to treat separately all the ^-intervals of the form, 
[l/(m + l),l/m], giving rise to an ever-growing set of wave number regimes but, in fact, 
it suffices for the results we wish to prove, to treat the last of these intervals, [l/m*,2A^], 
separately and combine all the other intervals into a single regime. 

So we define wave number regime III to correspond to the interval 1/m*] and wave 

number regime IV to correspond to the interval [l/m*,2A^]. These ranges are equivalent to 
•Sat < s < Sm* and Sm* < s < ir/h respectively. We will obtain a uniform bound on the 
magnitude of U^{s) for ^ G [l/(m + 1), 1/m] in regime III and an asymptotic value for the 
magnitude of U^{s) in the regime IV. 

Table [T] summarises the wave number regimes giving their boundaries in terms of both ^ 
and s. (Only the regimes for positive wave numbers are shown. By symmetry, the regimes 
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for the negative wave numbers are just the negatives of these intervals). The results of 



Wave number regimes 


Reg 


ime 




^max 




Smax 


I 

II 

III 

IV 




1 

N 
1 

m* 


1 

N 
1 

m* 

2A2 




SN 
Sm* 


SN 

Sm* 

n 
h 



Table 1. Limits of wave number regimes in terms of ^ and s. Here, r < 1/3, 
= 2A^ sin'^{h^^'^ /2), s„ are defined by 1/n = 2A^ sin^(sn/i/2) for n = and 
n = m* , and m* > 1 is the smallest n for which there exists such Sn 



carrying out asymptotic expansions in the four wave number regimes are summarised in 
the following Lemma which gives the behaviour of the error between the two transforms, 
E{s) = U^{s)-u{s,l)- 

Lemma 1. The following formulae give the errors in the Fourier transform in each wave 
number regime (see TableU^. 

(1) In wave number regime I, 

E{s) = u{s, 1) - ^A^.e + 1a2,4^ h\l + o(l)). 



(2) In wave number regime II, 

E{s) = oih'^) 

for any q > 0. 

(3) In wave number regime III, 

^ ((rn + l)!)2(iV-m-l)! 
' ^ ^' - {2m + 2){N + m)\ 

whenever ^ = 2A^ sin^(^) € [l/(m + 1), 1/m] and, everywhere in this regime, 

E{s) = 0(/i2™*+i) 

for some m* > 1. 

(4) Finally, in wave number regime IV, there exists a positive, continuous, and hence 
bounded, function S{-,m*) on [l/m*,2A^], such that 

\Eis)\=SiC,m*)h^'^l\l + 0{h))). 
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Proof. The key features of the proof are given here while further details are included in the 
Appendix. 

We first consider wave number regime I. Inspection of equation ^ shows that all the terms 
in the numerator are positive if ^ < < l/(-/V — 1). We can then write log(C/^(s)) as 

N-l 

log(C/^(s)) = J2 (log(l - ^0 - log(l + kO) - log(l + NO (4) 

k=l 

and expand the component terms, lo gfl =fc by f i rst ex panding ^ in powers of h. The 
calculations are very similar to those in lCarter &: GilesI (|2007l ) and the details of the expansions 
are left to the Appendix. We find that 

log{U^{s)) = -\s' + (^.^ - j^Xh' + h's^^ + Ois'h^). (5) 

we find we require that the second term converges to zero as /i — )■ and the last term is 
asymptotically dominated by the second term. These conditions are satisfied by requiring 
that s < h~'^ with q < and this condition on s defines the upper limit of the wave number 
regime I. 

1 2 

We then have, with u{s, 1) = e~2'^ , 

- u{s, 1) = u{s, 1) Q^s' - j^X's' + ^A^/) h\l + 0(1)) 
as we set out to prove. 

We now consider wave number regime II. 

The analysis is again similar to that in ICarter GilesI (|2007l ^. The lower limit of wave 
number regime II is 0{h~^) (which corresponds to the upper limit of wave number regime 
I) and the upper limit of the regime corresponds to the first zero of U^{s), i.e., s = 
{2/h) sin-i ^/{l/2X^{N - 1)), which is 0{h-2). In this regime each of the factors (1 — 
m^)/(l + (m + 1)^) in (s) is positive and decreasing as s increases. 

Consequently, the value of C/^(s) will be less than or equal to its value at s = h^^ for 
< r < 1/3. However, at s = h~'^ we have, for any q > 0, 

U^(s) ^€2^ / 1 X^\ ,ofT-'i^\ 1 

hi ~ W 

as /i — )• 0, so U^{s) = o{hi) at s = /i^3 and so we have U^{s) = o{hi) throughout this regime. 



Next we consider the high wave number regime which is split into two parts. 
It is in wave numbe r regimes III and IV that our analysis differs significantly from that in 



Carter &: GilesI ( 2007 ) where there are repeated factors in the corresponding expression for 



{s), which results in a simpler expansion in terms of 1/^. In our case, U^{s) is a product 
of many more factors, so the analysis and expansion is more complicated. However, we can 
obtain our key results by breaking the high wave number regime into two parts and we will 
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find that we only need to perform an expansion in regime IV. 



For wave number regime III, which extends from sn = (2//i)sin ^ {1/2\^N) to Srn*) 
corresponding to ^ running from 1/N to 1/m*, with m* > 1, we investigate the behaviour of 
U^{s) for values of s corresponding to values of ^ in the intervals [l/(m + 1), 1/m], where 
m* <m<N-l. We write U^{s) as 

f l(l-0(l-2e)...(l-mO 1 f (l-(m + l)0...(l-(Ar-l)e) ) 
^' \(l + ^)(l + 2e)...(l + (m + l)e)j 'l (l + (m + 2)e)...(l + iV0 /• ^' 

For the range of values of ^ G [l/{rn + 1), 1/m] being considered, the left-hand expression is 
a product of non-ncgativc factors, (1 — z/^)/(l + (u + 1)^), for u = 0, ...,m,, each of which is 
decreasing in ^, so the left-hand expression is bounded by the value of the expression evalu- 
ated at ^ = l/(m + 1), i.e., ((m + l)!)2/(2m + 2)! 

Similarly, the right-hand expression is a product of factors (1 — + {i^ + 1)0) fo^' 

u = m + 1, ...,N — 1, which are all negative and whose magnitude is increasing with ^. So 
the magnitude of the right-hand expression is bounded by the magnitude of the expression 
evaluated at ^ = 1/m, i.e. {N -m - l)!(2m -I- 1)!/(A^ + m)! 

Therefore, the magnitude of U^{s) for l/(m -|- 1) < ^ < 1/m is bounded by 

_ ((m + l)!)^(iV-m-l)! 
^^•"^ - (2m + 2)(iV + m)! ' 
and we consider the behaviour of this expression as m decreases from A" — 1 to m* . We have 

WN,m-i _ N^-m" 
WN,m rn{m + 1) ' 
and this expression is less than one for 

N'^ -m^ < m{m + 1), 

which corresponds approximately to m > N/^/2. This implies that, as m decreases from 
N — 1, Wiv,m first decreases, reaching a minimum in the vicinity of N/\/2. So the maximum 
of |L'"^(s)| will occur either at m = TV — 1 or at m = m*. 

We now compare the values of WN,m at these values of m, by calculating the ratio 

WN,m* _ {{m* + [N -m* -1)\{2N)\ 

Wn,n-i (2m* + 2) {N + m*)\ {N\f' 



By Stirling's formula, we have the approximation 

{2N)\ (2Af)^^+^e-^^ 2^^+^ 

(7V!)2 ~ ^^^^N+\^-NY ~ JV5 
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and we see that the exponential term dominates the other powers of N and WN,m* /Wn^n-i 
oo as — )• oo. So, for large enough N, W^r^j is maximised at m = m* and it then follows 
that the expression 

(2m* +2) {N + m*)\ " ^ ) - ^^'^ ) 

gives a uniform bound on |C/^(s)| in wave number regime III. 

Finally, we consider wave number regime IV. 

In this regime, we will find that we must perform an expansion of U^{s) from 
we write 



where 



so 



\U^is)\ =M{C,m*) 



(l-(m* + l)0...(l-(iV-l)0 



(1 + + 1)^) . . . (1 + (AT _ 1)^) 

(i-e)(i-2e)...(i-m*e) 



(i + e)(i + 2e)...(i + m*)o 



and 



logf|[/^fs) 



log(M(^,m*))+log 



[l 



+ log((l + iVO"') (8) 



+ 5(»n*+l) )•••(! + £,{N-l) ) 

and we seek an expansion of the left-hand side of this equation in terms of 1/^ to find, as 
described in the Appendix, that 



\U^{s)\ = 5(e,m*)/i(2/5+i)(i + o(/i)), 



(9) 



where 5" is a function of £^ and m*, which is positive and continuous on [l/m*, 2A'^], as we 
wanted to show. 

□ 

4. Error contributions from the various wave number regimes 
The error in the numerical scheme at Xj is found by using the inverse Fourier transform 

1 



Ej = —— / E{s) exp{—isxj) ds, 
2vr Js=-^ 

h 

where E{s) is the error in the Fourier transform. So decomposing the error, applying the 
inverse Fourier transform and using symmetry, we can write 
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where 



E) = - I E{s) cos(sxj) (Is, (10) 



vr 



el. 



where /j is wave number region for i G {1, 1 1, 1 1 1, IV}. We then have the following Lemma 
which makes use of the results from Lemma [1] in Section 3. 

Lemma 2. With E'^ = {Ej)jQZ defined in ilO\) . the contribution to the error of the numerical 
scheme from the four wave number regimes is given by 

E^ = 0{h^), 

E^^ = o{h'^) for any (7 > 0, 

E^^^ = Oih^""') for some m* > I, 

?IV _ n fe^ 1 



E'^ = O 



Proof. Again, only key features of the proof are given here while further details are included 
in the Appendix. 

In wave number regime I, the analysis is similar to that in ICarter Gilesl (|2007h . The key 
observation is that the term n(s, — ^A^s^ + |A^s^) has a readily identifiable inverse 

Fourier transform, and we deduce that the error contribution is 0{h'^). 



Li wave number regime II, we again follow the analysis in (jCarter &: Gilesl . 120071 ) to con- 
clude that the error contribution from this regime will be dominated by h'^ for any g > 0. 

In wave number regime III, we can write 

\U^{s)cos{sxj)\ds < \U^{s)\ds, (11) 

-N-l •'0 

which is TT/h ■ 0(/i2"^*+i), i.e., ©(/i^"^*). 

Finally, we consider the error contribution from wave number regime IV, i.e., where Sm* < 
s < Tr/h. 

We need to determine the behaviour, as — )• 0, of 

S{^,m*)h^^~^^^ cos{sxj) ds, 

the magnitude of which is maximised at j = (this is because S is positive in this wave 
number regime). So we only need to consider the value of 

F{h) = / " 5(e,m*)/i^^+t) ds. 



From ^ = sin2(^), we find that 
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SiC,m*)hl 

2 

As the lowest (i.e., dominant) order in h from /i? occurs when ^ = 2X^ we expect that this 

1 

integral will be close to 0{h~>^) and using careful, but elementary, analysis (see the Appendix) 
we find that 




hJ^ V Vlog l/h I 



and so the error contribution in regime IV is found to be O (h^ / \ log ^ ) . 



□ 



As a consequence of Lemma [21 we obtain Theorem [H as stated in the Introduction. This 
result follows from the observation that when A < l/\/2, we have m* > 1 and 1/A^ > 2, so 
the convergence behaviour is eventually dominated by wave number regime I and the errors 
(in the maximum norm) will be O^h?). However, if A > l/\/2) then m* = 1 and 1/A^ < 2, 
so the numerical error of the scheme is dominated by wave number regime IV. The scheme 
still converges but the errors will now be 0{h^^^ /-y/log(l//i)). The logarithmic factor in this 
expression is slow growing and the error will be close to 0{h^^'^^) which is worse than 0{h?). 

So our result is that we always get convergence in the transformed numerical scheme 
whereas, without the time change, the errors eventually increase as /i — >• 0, with A controlling 
how small h has to be for the divergence to manifest itself in practice. 



5. Numerical results and complexity considerations 

Numerical experiments have been performed to investigate the behaviour of the time- 
changed scheme used to solve the heat equation with delta-function initial conditions. The 
analysis used T = 1 and the space domain was truncated at x = ±10. The numerical scheme 
was run with a largest value of A; = 0.01, i.e., with 100 time steps in the transformed coordi- 
nates. The grid was repeatedly refined by dividing both h and A; by 2 and the calculations 
were stopped when the number of time steps reached 3200. 

The error (in the maximum norm) of the scheme was analysed by comparing the numerical 
results with the exact solution and the slope of the log-log plot of the error against h was 
calculated. The slope was then plotted against the value of A in Figure 1, along with the 
h exponent for the error derived from the analysis above, where we ignore the effect of the 
log(yl//i) term on the behaviour of the error in wave number regime IV and just use the 
expression min(2, 1/A^). 

We see a good match between theory and experiment with the larger divergences being 
seen where A is closest to the critical value of l/\/2. Figures 2 and 3 show the log- log plots 
of the errors versus h for the cases A = 0.5 and A = 1.0, either side of l/\/2. 

It is not possible to determine experimentally the h exponent with complete accuracy for 
two reasons. 
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Figure 1 . Convergence order for the time-changed scheme as a function of A 
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Figure 2. Error plot for the time-changed scheme with A = 0.5 

Slope = -2 




Firstly, as /i — )• 0, the time required for each solution increases rapidly so, for compari- 
son purposes, all the calculations were terminated when the number of time steps reached 
3200. As the overall error of the scheme is a mixture of errors with different h exponents, 
the asymptotically dominating h exponents will not be reached because the calculations are 
terminated early. In particular, for A slightly above l/\/2, the contribution to the error from 
the highest wave number regime will be limited and lower values of h need to be reached 
before this error clearly dominates the other errors. Note that our analysis does not identify 
the weighting factors for the errors of different order in h. 

Secondly, the additional term, ^J\og{l/h), in the error behaviour in wave number regime 
IV will distort the error plot and tend to increase the apparent h exponent to some extent. 
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Figure 3. Error plot for the time-changed scheme with A = 1.0 

Slope = -1.05 




1/h 



We have also carried out a comparison of the error performance between the time change 
scheme and the Rannacher scheme. For A < 1 /\/2, both schemes will show quadratic conver- 
gence and, for small values of h, we can write 



U^^ {s) - u{s, 1) = n(., 1) ( j^s' + ^A^.^ - ^Xh' )hHl + 0(1)) 



and 



U^cis) - n{s, 1) = n(., 1) { j^s' + ^A^.^ - ^A^.^ ) h\l + o(l)), 



where (s) and U^(j{s) are the Fourier transforms for the Rannacher and time-changed 
scheme in wave number regime I. 

For m > 1, the mth derivative of the cumulative Normal distribution has the Fourier 
transform given by 

iVM(s) = (is)(""-^)e-^'/2 
so that the Fourier inverse of (fs)(™~^)e^'^ is —^N^^'^ix). We can then estimate the ratio 

V 27r 

of the errors of the two schemes (at x = 0) as 

Er 



3 + 1^^) 



15g^A^ 



Etc 3(^ + iA2)-15^A2 
by evaluating iV(^)(0) = 3 and iVC^)(0) = -15 from 

1 



and 



iV(7)( 



15x^ + 45x2 - 15)e-^ 
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The ratio Er/Etc has been plotted as a function of A in Figure [Hand the results compared 
to the empirical observations. So for values of A close to the critical value, there is a modest 
reduction in the error of the numerical scheme due to the time change method. 

Note that for A > l/\/2 the errors will no longer be comparable as the order of convergence 
of the time-changed scheme is no longer quadratic. 



1.5 



Figure 4. Plot of ratio of Rannacher to time change error 
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We are also interested in the computational efficiency of the scheme, so we have investigated 
the relative performance of the two schemes when both are subject to a constraint on the 
computational cost, C. We can write C ~ c[\/k){\/h) for some constant c and we then can 
approximate the errors as 

Er = 3(1/24 + AV8)/i2 - (15AV96)/i^ 
= (l/8)/i2 + (21/96)/t2 



and similarly 



Then 



and similarly 



Etc = (l/8)/i2 + (3/96)A:2. 

Er = hk{{l/8){h/k) + {2l/96){k/h)) 
= (c/C)((l/8)(l/A) + (21/96)A) 



Etc = (c/C)((l/8)(l/A) + (3/48)A). 
The errors Er and Etc will be separately minimised as functions of A at 



X*R = x/(l/8)(96/21) ~ 0.756 



and 



^/(iTsKWsy = \/2. 
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Figure 5. Plot of ratio of Rannacher to time change optimal error at the same cost 

Rannacher v time change errors for fixed computationai cost 



nJ 0,25 - 

E 



2 c 



Note that for the time change the analysis is only valid for A < l/-v/2i so the minimum 
occurs at \|^/2. 

The plot (Figure [5|) of the errors for A < l/v^, with the fixed computational cost, shows 
that the time change scheme has the lowest achievable error for the value of A = l/\/2 and 
this error is, in turn, better than the error of Rannacher scheme at the value Ao. 



6. The time change applied to the Black-Scholes equation 

In this section, we describe numerical results obtained by applying the time change to the 
Black-Scholes equation 

^+r''^'S+'^l^-'^=°' ^>0'*^(0'^]- (13) 

This is the equation that can be used to calculate, for example, the price, V ^ of a European 
option. Here, S is the price of the underlying, a is the volatility, r is the risk-free rate, and t 
is time. For a call option, the terminal condition is given by the payoff, i.e., max(5 — -fC, 0), 
where K is the strike of the option. Practically important sensitivities of the call price to 
the initial asset price, used in risk management, are given by the delta. A, and gamma, F, 
defined as 

"-'as 

and 

respectively. 
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Our investigation of convergence for the heat equation was based on Fourier analysis that 

relies on the PDE having constant coefficients whereas the Black-Scholcs PDE used to price 
a European option does not have constant coefficients and so our analysis is not directly 
applicable in this case. 

Nevertheless we have implemented the natural time change method to study the conver- 
gence behaviour of the gamma calculated with the modified numerical scheme. 

In this context the natural analogue to the solution of the heat equation with Dirac delta 
initial conditions would be the PDE satisfied by the gamma, 

f + + + 2.^)1 + (r + .')r = o, 

where the terminal condition for this equation is given by the second derivative of the payoff 
with respect to S, i.e., T{S,T) = 6{S — K). Here we have also assumed that a and r are 
constant. 

This PDE could be used to test the performance of the time change approach but a more 
useful approach (as it reflects industry practice) is to apply the time change directly to the 
Black-Scholes equation and then use finite differences to derive the values for the gamma from 
the calculated option values. (Note that the scheme was also tested with the gamma PDE 
and gave identical results). 

We first note that if we can write our PDE in the form 

■ 

where r = T — i is the time to maturity and where the spatial differential operator C has 
time-independent coefficients, then the transformed equation for t = ^Jr is simply 

dV 

~2tCV = 0. 

dt 

If the coefficients in C are time-dependent, then these need to be rewritten in terms of i. 
This simple change to the equation results in very minimal changes to the solver code. The 
calculation done for a single time step in the original variables can be written as 

(7 + L)y"+i = (7 - L)y", 

where L is a matrix dependent on k, h, S, a and r while is the solution vector at time 
step n. In the case of constant coeflacients, and in the time-changed variables, this calculation 
simply becomes 

(7 + 2 fn+iL)F"+i = (7 - 2 inm"", 

where serves as numerical approximation to F(-,t„) = V{-,T — tn^). In particular, the 
transformed Black-Scholes equation is 

with the payoff function as the initial condition. 

The transformed equation for V was solved using the Crank-Nicolson scheme and the 
gamma calculated by finite differences. The errors in gamma were computed from the analytic 
values and the convergence rate of the numerical error as a function of A was analysed. 
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Although the scheme does not have the constant coefficients required for the application 
of Fourier analysis, the errors do behave in a very similar way to those for the heat equation 
(see Figure 4). 



Figure 6. Convergence order for the gamma calculated from the time- 
changed Black-Scholes scheme function of A 
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The volatility, cr, was set to 20% and the strike, K, was set to 100. The risk-free rate, 
r, was set to 5%. The results of the analysis showed that there was a critical value of A 
below which the error in the maximum norm was 0{h?) and above which the h exponent 
declined in the same manner as for the heat equation. Also shown in Figure 4 is the plot of 
min(2, l/((T^i^^A^)) and it is clear that this gives a reasonable description of the behaviour 
of the h exponent although it is not quite as good as with the heat equation itself. 

We can justify the critical value of A by noting that, in the Black-Scholes equation, the 
factor plays the same role as the factor ^ in the heat equation and we also observe that 

if we had written the heat equation in the form ut = \D Uxx, for some constant diffusivity D, 
we would have found the critical value of A to be l/y/2D. 

Typically, the worst behaviour of a numerical scheme for the Black-Scholes equation with 
delta function initial conditions is seen at times close to maturity and in the vicinity of the 
strike and so the strike is an appropriate parameter to use to represent the coefficient of 
which then determines the critical value of A. 

Even without the precise results available from Fourier analysis, we can still attempt a 
partial explanation of how the time change improves the convergence of the Crank-Nicolson 
numerical scheme. One consideration is that the truncation error of the numerical scheme 
has a leading order term proportional to the third derivative of the European option price 
and the behaviour of this derivative is much improv ed by t he tim e change. 

The theta, 0c, of the call option is given by (see Hand ( 2007 ). p. 64) 



where C is the solution to (fT3|l with C{S,T) = max(S' — K,0), N{x) is the cumulative 
Normal(0,l) distribution, n{x) is its first derivative (the Gaussian function) and 
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{\og{S/K) + {r±\a^){T-t)) 

Clearly Qq will become singular at the strike as t — t- T and any higher derivatives will also 
be singular there. 

However, if we change the time variable to t = \/T — t, we find that 



dt dt \ 2t ^ ' 

= -Sn{di)a + 2irKe-''^^N{d2) 
and @c is then not singular at t = 0. 

Carrying on to the second 'time' derivative we see that 

^ = -Sn'idi)a^ + 2rKe-''\l - 2 t r)iV(d2) - 2 t rKe-''\{d2)^, 
o¥ at dt 

where we also have 

5di,2 _ log f ^ r ± 4 



dt af^ a 

so there are potentially singular terms to consider, i.e., those featuring inverse powers of t. 

However, S ^ K then both di and d2 — ?• ±oo as t — )• so that the exponential factors in 
n{di) = (l/\/27r) exp {—df/2) will eventually overwhelm the inverse powers of i. 

Finally, ii S = K the term log(S'/i^) vanishes and the derivative is non-singular at the 
strike. 

It is clear that this analysis will extend to the higher derivatives of C because the negative 
powers of t will now only arise from differentiating N and will then only appear as multipliers 
of the Gaussian function. These powers will only be present for S ^ K where the Gaussian 
function will decay rapidly and overwhelm the inverse powers. So the higher derivatives will 
not become singular as t — )• 0. It is also clear that the factor t in ([H]) causes the singular spatial 
truncation error, which contains the fourth derivative of V with respect to S, to become more 
well-behaved for small t in the vicinity of the strike. So the time change has, to some extent, 
regularised the behaviour of the truncation error. 

Although it is not entirely clear how to deduce convergence properties of the transformed 
numerical scheme directly from this, and particularly convergence of derivatives of the solu- 
tion, it seems plausible that the improved regularity of the truncation error is advantageous 
when using initial conditions that are not smooth. 

7. American options 



Forsyth & Vetzal (|200i) investigate the convergence properties of Crank-Nicolson time- 



stepping when used to calculate the price of an American put by a penalty method. The 
penalised Black-Scholes equation is given by 

dV 1 d^V dV 

'dt + 2'"''^' + ^^55 - + Pmax(max(if - S,0) - V, 0) = 0, 

where p > is a penalty parameter. This non-linear equation is solved with the terminal 
condition set equal to the option payoff. The penalty term comes into play when the solution 
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violates the requirement that V should not fall below the payoff. As the penalty parameter 
p increases, the solution of the equation converges to the value for the American option 
satisfying 

^ dV 1 d'^V dV \ 
max f — + 2^'^'^ + rS-^ - rV, maxiK -S,0)-V]= 0. 

For the present purposes, we are interested solely in the convergence of the Crank-Nicolson 
scheme and choose p large enough to make the penalty solution numerically indistinguishable 
from its limit. 



Forsyth & Vetzal ( 20021 ) implemented their solver using the Crank-Nicolson scheme with 



Rannacher start up and demonstrated that as the time step was reduced with the ratio of 
time step to space step held constant, the error in the calculated at the money (ATM) value 

3 

of the American put tended to 0(/i2) behaviour. Their analysis of the problem suggests 
that non-uniform time-stepping might restore second order convergence so they went on to 
describe and implement an adaptive time-stepping procedure that did indeed yield second 
order convergence for the option value. 

Their results suggest that the simple time change method described above might also be 
able to restore second order convergence. 



Ratios of successive differences 



A = 0.0125 



M 


N 


Value 


Delta 


Gamma 


3200 


640 


3.98 


3.92 


3.91 


6400 


1280 


3.97 


3.97 


3.94 


12800 


2560 


3.97 


3.97 


3.96 


25600 


5120 


3.99 


3.97 


3.99 



Ratios of successive differences 
A = 0.0250 



M 


N 


Value 


Delta 


Gamma 


3200 


320 


3.84 


4.01 


3.80 


6400 


640 


3.99 


3.93 


3.96 


12800 


1280 


3.94 


3.93 


3.89 


25600 


2560 


3.96 


3.97 


3.98 



Ratios of successive differences 
A = 0.0500 



M 


N 


Value 


Delta 


Gamma 


3200 


160 


3.85 


3.65 


2.07 


6400 


320 


3.78 


3.87 


2.08 


12800 


640 


3.89 


3.87 


2.09 


25600 


1280 


3.89 


3.90 


2.09 



Table 2. Ratios of successive differences between numerical solutions for suc- 
cessive refinements for V, A and T; M the number of steps in S, N the number 
of steps in t; 0- = 0.20, r = 0.05, p = 10^ 
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The penalty code was modified to incorporate the time change and numerical results were 
generated. Second order convergence was observed for the option value and delta for the 
ATM option and it was also observed that, for sufficiently small A, the gamma showed second 
order convergence. Table [2] shows the ratios of successive differences for mesh refinements 
with A equal to 0.0125, 0.025 and 0.05. As with the European option discussed above, and 
with the strike, K = 100, and the volatility, a = 20%, the critical value of A appears to be in 
the vicinity of l/{aKV2) = 0.035. 

It was also observed that the plots of the calculated gamma revealed i nstab ility around 
the exercise boundary, behaviour which was also seen in iForsvth fc Vetzall (l2002h . where the 
authors found that their adaptive time-stepping approach resulted in instability of the gamma 
at the exercise boundary unless fully-implicit time-stepping was used. This behaviour is a 
separate phenomenon independent of the short-time asymptotics addressed here. 

In the previous section we suggested that the improved convergence properties of the Crank- 
Nicolson scheme might be related to improved regularity of the third time derivative of the 
value function of the European put in the vicinity of the strike and close to maturity. 

We would like to be able to perform a similar analysis for the American put but, of course, 
nosimple closed form expression for the value function is known. However, Mallier fc: Alobaidil 
(|2004l ) do present an asymptotic expansion for the value function which can be decomposed 
into the value function of the European put plus the sum of an asymptotic series representing 
the e arly exercise premiu i n. 

In iMallier fc Alobaidil (j2004[ ) , some changes of variable are used to define the final form 
of the expansion. We write S = Ke^ and ^ = x/{2^JtI), where Tg is the rescaled time to 
maturity, equal to 2(T — t)/a'^. The transformed value function is then given by 

where P{S, t) is the value of the American put as a function of the original variables. We can 
then write that 



v{x,Ts) = v^{x,T,) + r;/2(-logr.)— Fr(0, (15) 



=2m=l 

where are functions of the similarity variable ^ and v^{x, Ts) is the value of the European 
put e xpressed in the new v ariab les. 

In iMallier Alobaidil (j2004l ). only a few of the funct ons are fully identified and 



discussion of those is restricted to m = 1. Those with n > and m = 1 are, in principle, 
tractable but only partial results are known for the cases where m > 1. This limits our ability 
to draw conclusions from the behaviour of the derivatives of the early exercise premium. 
Nevertheless, we will illustrate some aspects of the behaviour of the derivatives by describing 
a single term of the early exercise premium and its form and behaviour after the time change. 
For example we find that the leading order coefficient of the early exercise premium is given 

by 

F^iO = tx{ie^v{-em/V^+{2e + l)erfc(C)) 

for some constant fi. 

After the time change this early exercise term in (fT5|) is 

G = ^,?{-2\ogi)-^Fl{l^), 
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where rj = x/{2t). At the strike, we have 

n{2i{-2logiy^ - 2i{-2logiy'^), 

/x(2(-2 log i)-^ + 4(-2 log i)-^) + 8(-2 log iy^) 
and 

-J- = ^(4(-21ogt)-' + 16(-21ogt)-3) + 48(-21ogt)-4). 

Acknowledging that this calculation only addresses a single term of the expansion of the 
early exercise premium we see that there is an improvement in the regularity of the derivatives 
of the put but unfortunately this does not obviously extend to the third derivative terms such 
as d^G/dF. 

So, although we are encouraged by the empirical performance of the time change method 
applied to the American put, further analysis will be required to explain fully the improvement 
in convergence. 



dG 
'dt 



8. Conclusions 

In this paper we have presented the analysis of a simple and natural change of time variable 
that improves the convergence behaviour of the Crank-Nicolson scheme. When applied to the 
solution of the heat equation with Dirac delta initial conditions, the numerical solution of the 
time changed PDE is always convergent with rate of convergence determined by the ratio. A, 
of time step to space step (held constant during grid refinement). 

This behaviour contrasts strongly to that of the Crank-Nicolson scheme when applied to 
the original PDE, which is always divergent and where A controls how small the time step 
must be before the divergence appears. A proof of the behaviour of the time-changed scheme 
is given and numerical results presented to support the theoretical analysis. 

Although the Fourier analysis used to prove this result cannot be applied directly to the 
Black-Scholes equation as it does not have constant coefficients, numerical experiments for 
the price and Greeks of a European call indicate that the time change method also leads to a 
convergent Crank-Nicolson scheme for this problem, without the need to introduce Rannacher 
steps, and quadratic convergence can be obtained if A is chosen appropriately. 

The change of convergence order at the value of A which gives optim al complexity fo r given 
error tolerance is in line with the sensitivity of the error to A o bserved in iGilesfc Cartel (I2OO6I 



The present analysis does not support the conjecture made in iGiles Cartel! (|2006l ^. however 



that the convergence of the Rannacher scheme with non-uniform (square-root) time-stepping 
is never of second order, the argument being that for small initial time-steps the smoothing 
effect diminishes. In fact, we see here that under square-root time smoothing is not necessary 
at all. 

The experimentally observed second order convergence for the time-changed scheme ex- 
tends to the computation of American option values which satisfy a linear complementarity 
problem with sin g ular b ehaviour of the free boundary. These results are in line with those of 
Forsyth &: Vetzal ( 20021 ) for a time-adaptive Crank-Nicolson scheme with Rannacher start-up. 



It can be seen as an advantage of the time-changed scheme proposed here that it provides 
a remedy for both the reduced convergence order of the Crank-Nicolson scheme for European 
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and American option values and for the divergence of their sensitivities by a single, simple 
modification. 

Future work in this area will investigate other applications of the time variable change to 
option pricing and will consider whether the method offers efficiency gains over the standard 
approaches. 
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Appendix A. Proof of Lemma [H ([5]) and ([9]) 

We analyse the behaviour of the Fourier transform U^{s) in each of the four wavenumber 
regimes. 

Proof of dl]): We start from (gD, 

N-l 

log{U^{s)) = (log(l - k^) - log(l + kO) - log(l + m). 

k=l 

We seek to expand the logarithmic terms. We begin by first expanding ^ = 2A^sin^(^) in 
powers of h to obtain 

l-k^ = l- -^h" + -^h^ - + -^^^k = + 

i=0 

where \9k\ < 1- Setting 6 = h'^, we define 

3 

= 1 - A;^ = ^ Ai6' + ^4^^ 

i=0 

and 

/(<^) = iog(i-A:e) = iog(5(5)) 

and we can then expand f{5) in terms of S to obtain 

and we have to ensure that the remainder Z/^ is well-behaved. We do this by analysing the 
behaviour of 

where 

p{5) = g{5fg^*H5) - igi6)g' {6)g(-^H^) + 6{g'i6))^gi6)g"{6) - 3ig"i6)fgi5) - 6ig'i6))\ 

Then g{6) — )• 1 as /i — )• if we ensure that s < h^"^, where m < ^, because, for example, 

, , kXh^h^ NX^s^h^ Xs^h 

\Ai6\ = < = > 0. 

' ' 2 2 2 

Then h can be chosen small enough (or N large enough) so that 

\1 + Ai5 + A26^ + Asd^ + Ai6^\ > ^ 

so the denominator of /^'^^ (6) is bounded away from zero. Next we can show that 

\g^'\6)\ < ai\Ai\ 
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for some constants, at for i = 0, 1, 2, 3 and 4, because of the constraint on s. We can finally 
conclude that 

|/(^)(#)|(5^ < (2a4|^4| + 32aia3|Ai||A3| + 48a?a2|Ai| V2I + 2Aal\A2f + 96af\Ai\^)h^ 

and it follows that we can estimate asymptotically the remainder by just retaining the term 
proportional to lAil^/i*^ which is, in turn, proportional to k'^X^s^h^. 

By changing the sign of k, we can obtain a similar expression for log(l + k£^) hence find 
that 

log(l - kO - log(l + kO = -kXh^h^ + —kX^s'^h^ - —kX'^s^h^ - —k^Xh^h^ + Wkk^s^h^ 
B\ Bv ^2 360 12 

where \Wk\ < K for some constant, K, and so 

/AW 1 ;^2^4^4 _ J_xh%A N + -X^s^h^N^ - ^X^S%'N' 
V 2 24 720 J 8 48 

- ^X^s%^N^ + 0{s^h^), 
because (EtZi k^)^^s^h^ = 0{N''Xh%^) = 0{Xh^h^). So 

log(f/^(.)) = ^-XW + ^^XW - ^XW^ (i(iV - l)iV + ^) 

as the remaining terms are dominated by 0{s^h^). Substituting N = we get 

log(&"(»)) = -lo^ + (ij.' - + Ia^,') - ( + 1a^,» + 1a^.») + 0(»»tf ) 

and hence 

because the s^h'^ term of the expansion is dominated by the term 0{s^h^). 

Proof of do]): We start from ([HD, 
log {\U^{s)\) = log(M(e, m*)) + log( i|"^|"""||"^U + log((l + iVO-^), 

[ + C(»n* + l)J • • • U + ^(N-1)) ) 
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and write 



k=N-l 



1 \ , / 1 



We can expand the logarithmic terms for ^ G [l/m*,2A ] to get 



k=N-l , , , ^ k=N-\ / oo 



E iog(i-j!,)-.og(i + i U E EM-E 



(-1) ^ (-1?+^ 



fe=m*+l ^ ' fe=m* + l \ j=l ^ ^' j 

" ~ ^ (2i + l)^2i+i I H pi+ij 

j=0 ^ \A:=m*+l / 

and, for any N , this series converges absolutely in 1/^ as we have just rearranged a finite sum 
of absolutely convergent series. We can write this as 

k=N-l / , ^ 2 ( l\ °° 2 / k=N-l ^ \ 



E =-| E r -EottW E 



fe=m*+l 3 ^ a \fc=m*+l / j=l \k=m*+l j 

^ - +A(C,m*) + B(C,iV), 



\fc=m*+l 



where 



and 

oo 2 / °° 1 \ 

and we must establish the convergence and size of the functions A(^,m*) and B{^,N) (see 
below). 

We also have 



■ iog(i +N0 = - log m - log (1 + = - log - f; 

= -logiV^-^ + C(C,iV), 
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where 
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,=2 ■^■(^^)' 



SO that 



' k=N-l 



log(|C/^(^)|) =log(M(e,m*))-| ^ i j +A(e,m*)+i?(^,iV)-logiVe-^+C(e,iV) 

2 ( \\ 1 

= log(M(e,m*))-- J] _ +A(e,m*)+i?(C,iV)-logiVe + ^+C(e,iV) 

\k=m*+l / ^ 

= P(e,m*)- - ( ^ -1 + D(e,iV)- log iV + 0(iV-i), 



where 



^k=m* 



and 

Then we can write 



P(e,m*) = log(Af(e,m*))-||^^ +A(e,m*)-loge 

i)(e,iv) = s(e,iv) + c(e,iv). 



log(|[/^(s)|) =P(^,m*) - I (logAf + 7 + O(iV-i)) +D(e,A^) - log 7V + 0(iV-^) 

= Qit m*) - (I + l) log ^ + 0(iV-') + ^(e, N), 
where 7 is the Euler-Macheroni constant (see Trij . I2OIII ) , and where 



Q(e,m*) = P(e,m*)-|7. 

Then using N = l/{hX) and assuming for now D{^,N) = 0{N~^), which we will show later, 
we get 

log = R{^,m*) + (I + 1) log/i + 0{h), 

where 

i?(C,m*) = Q(C,m*) + + 1) log A. 

Then we finally have 

=5(^,m*)/i(t+0(l + O(/i)) 

(16) 
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as we were required to prove. We note that, written out in full, we have 

-| (7+ E I 



and we note that S{^,m*) is continuous and hence bounded on [l/m*,2A^]. 

We still need to prove that A{^,m*) + B{(^,N) is a valid rearrangement and proceed by 
proving that A(^,m*) is convergent. It will then follow that B{^,N) is also convergent. In 
addition we will show that D{^, N) = B{C, N) + C{C, N) = 0{N-^). 



By a result in lTimoftd (|2005l ). we have, for r > 1, 

A:=oo 



V ^— = - 

^ /fc^'^+i 2r{m* + Om^Y'' 



for some 9f)-i* ^ . Therefore 



k=m*-\-l 



i^(c,m*)i < x; 



^ (2j + 1) e^J+i 2j{m* + i; 



2j 



1 



00 

2i 



where = l/(^(m* + ^)) < 1 and so A{£^,m*) is convergent. It follows that B{^,N) is con- 
vergent because A{(^,m*) + B{£^,N) is convergent. 

Next we find out how B{S^,N) depends on A^. We have 

°° 1 1 1 
\BiC,N)\ < 2g(2i + 1)^2^+1 2j(iV- 1 + 1)2.- 

^ 1 1 N^^ 

^4^(2j + l)(iV^)2.+i2j(A^-i)2^' 

and as N/ [N — |) < 2, we have 

00 2j+l 1 °° / 2 \ 

m.N)\ < 2ivX: + 1) (iv;).,.. — < I^'^S (ivi ' 



2j + l 
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Then because we have 2/^ < 2m* , we find that for sufficiently large A^, we have 2/N£^ < c < 1 
for a constant c and the geometric series converges and we find that 

and similarly for C. It then follows that D{^,N) = 0{N-^). 

Appendix B. Proof of Lemma [2], (ri2|) 

We want to show that 



d^ = 

where S{-,m*) is continuous. We make the substitution = \ — ^/2A^ to obtain 

I = 4X^ S*{z,m*)hy^^) dz, 

Jz=0 

where B = 1 — and where S*{z,m*) = S{^, m*)/\/l — is a continuous function on 
[Q,VB] as B <1. 

We want to show that this integral is concentrated around z = for small h and write 

I = h + h. 



where 



and 



I^ = A\^ I S*{z,m*)hy^^) dz 



2 = 



I2 = 4X^ S*{z,m*)hy^^-'') dz, 

Jz=y/A 

where A is chosen so that, asymptotically, Ii dominates l2- 

Firstly, we consider I2. For /i < 1, the second factor of the integrand is decreasing in z, so 

1 A 



I/2I <4A25;;,,/iVi^i-^ 

where is a bound for 5* on [^m*>2A^]. We now take 

A 1 



and note that 



( 1 I A 1 1/1 

log/.i^-^=-^^log- = -Wlog- 



so 



\l2\<4X^S*e I^V^^O as /i ^ 0. 
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Now we want to show that 



lim/i\/log| = lim |4A^/log^ I S*iz,m*)h^'^^) dz\ = Oil) 
h^o V h h^o \ \ h J^^Q I 



and we note that 

A = ^ ^0 as /i ^ 0. 

1 + ^b^TTTI 

We make a final substitution 



to get 



lim/iy;^ = 4A^ lim r ' 5*(Ar;/Vl^, m*) exp {- Y^y^^fJ^^^^ 



'??=o 

where the upper limit 



tends to infinity as /i — > 0. As >S'*(-, m*) is bounded and continuous at we see that 

/•'?='?* , / \ j"n=co 

SO we have /i -^log 1/h = 0(1) as /i — >■ 0. 



From this we see that Ii = 0(1/ <J log ^) and we also have 





lim 



= limA/log^e"^^/^ = 

/i-5>o V h 



so we have established that /i asymptotically dominates I2 and so / = 0(1/ ^/log 
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